Untitled

graficador <- function(a, b, l, f) {
    x <- seq(a, b, length.out = l)
    y <- x
    z <- outer(x, y, f)
    z[z < 0] = 0
    op <- par(bg = "white")
    persp(x, y, z, theta = 30, phi = 30, expand = 0.5,
        col = "lightblue", ltheta = 120, shade = 0.75,
        zlab = f, xlab = formalArgs(args(f))[1],
        ylab = formalArgs(args(f))[2]) -> res
}

tabla <- function(dataframe) 
  # Genera tablas con código de colores para un dataframe
  {
    gt(dataframe, rownames_to_stub = TRUE) %>%
        gt_color_rows(columns = dimnames.data.frame(dataframe)[[2]],
            palette = "ggsci::blue_material", domain = range(dataframe)) %>%
        gt_theme_nytimes()
}

FUNCION DE PRECIO función relu

price_function <- function(p_1, p_2, k = 400)
{
    return(1 - (exp(p_1/k)/(exp(p_1/k) + exp(p_2/k))))
}
graficador(3, 1000, 60, "price_function")

FUNCION DE DECISON

dicotomic_function_ps <- function(c_ps, c_xb, l_ps = 0.8) {
  return((1 - c_ps) * (1 - c_xb) * 0.75 +
           c_ps * (1 - c_xb) * 0.9 +
           c_xb * (1 - c_ps) * 0.3 +
           c_ps * c_xb * l_ps)
}

dicotomic_function_xb <- function(c_ps, c_xb, l_xb = 0.9) {
  return((1 - c_ps) * (1 - c_xb) * 0.7 +
           c_ps * (1 - c_xb) * 0.4 +
           c_xb * (1 - c_ps) * 0.8 +
           c_ps * c_xb * l_xb)
}

CONTROL DE PRECIO

price_control_ps <- function(p, m = 2000, c_ps, c_xb,
    l_ps = 0.8) {
    m = m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps)
    return(0 + (p <= m) * (p>=0) *((-(p - m)/m)^2))
}
price_control_xb <- function(p, m = 2000, c_ps, c_xb,
    l_xb = 0.9) {
    m = m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb)
    return(0 + (p <= m) * (p>=0) * (((p - m)/m)^2))
}

FUNCIONES DE VENTAS

ventas_ps <- function(p_ps, p_xb, c_ps = 1, c_xb = 0,
    k = 400, m = 2000, l_ps = 0.8)
    {
    return(5e+07 * dicotomic_function_ps(c_ps, c_xb,
        l_ps = l_ps) * price_function(p_ps, p_xb, k) *
        price_control_ps(p_ps, m, c_ps, c_xb, l_ps = l_ps))
}
ventas_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_xb = 0.9)
    {
    return(5e+07 * dicotomic_function_xb(c_ps, c_xb,
        l_xb = l_xb) * price_function(p_xb, p_ps, k) *
        price_control_xb(p_xb, m, c_ps, c_xb, l_xb = l_xb))
}

FUNCIONES DE BENEFICIOS

beneficios_ps <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_ps = 0.8)
    {
    ventas_ps(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
        m = m, l_ps = l_ps) * (p_ps - 100)
}

beneficios_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_xb = 0.9)
    {
    ventas_xb(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
        m = m, l_xb = l_xb) * (p_xb - 100)
}

graficador(1, 1000, 60, "beneficios_ps")

graficador(1, 1000, 60, "beneficios_xb")

FUNCIONES DE RESOLUCIÓN

# la siguiente función calcula el equilibrio de Nash
# del problema de Bertrand
resolucion  <- function(c_ps = 0,
                       c_xb = 0,
                       k = 400,
                       m = 2000,
                       x0 = 0,
                       y0 = 0,
                       N = 100,
                       l_ps = 0.8,
                       l_xb = 0.9,
                       eps=10^-4)
  # x0 e y0 corresponde al inicio
  # N corresponde al número de iteraciones
  # usaremos por defecto N = 100 por funcionar bien
{
  # generamos una matrix que contenga todas las iteraciones
  z = matrix(0, nrow = 2, ncol = N + 1)
  z[, 1] <- c(x0, y0)# almacenamos el punto inicial
  i = 0
  l=TRUE#criterio de parada por diverger
  aux2=TRUE #criterio de parada al converger
  # iniciamos el bucle
  while (all(i < N, aux2,l))
  {
    i = i + 1
    #almacenamos el óptimo para p_p
    aux = optim(z[1, i], function(x)
    {
      -beneficios_ps(x, z[2, i], c_ps,
                     c_xb, k, m, l_ps = l_ps)
    })$par
    #comprobamos que no sea demasiado alto
    if (aux <= (m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps))) {
      z[1, i + 1] <- aux
    } else{
      l=FALSE
    }
    
    #almacenamos el óptimo para p_x
    
    aux = optim(z[2, i], function(y)
    {
      -beneficios_xb(z[1, i + 1], y,
                     c_ps, c_xb, k, m, l_xb = l_xb)
    })$par
    if (aux <= (m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb))) {
      z[2, i + 1] <- aux
    } else{
      l=FALSE
    }
    aux2=dist(z[1, i + 1]-z[1, i])>eps
  }

  #nos devuelve los beneficios de cada empresa
  # y los precios de las consolas del equilibrio de Nash
  return(list(
    'beneficios' = c(
      beneficios_ps(
        p_ps = z[1,
                 i + 1],
        p_xb = z[2, i + 1],
        c_ps = c_ps,
        c_xb = c_xb,
        k,
        m,
        l_ps = l_ps
      ),
      beneficios_xb(
        p_ps = z[1, i + 1],
        p_xb = z[2, i + 1],
        c_ps = c_ps,
        c_xb = c_xb,
        k,
        m,
        l_xb = l_xb
      )
    ),
    'precios' = z[, i + 1],
    'converge' = l
  ))
}

SOLUCIÓN

# La siguiente función calcula la p y la q 
# de un problema 2x2
solucion2 <- function(e1,e2) {
# importamos las funciones de python
np <- import("numpy", convert = FALSE)
nash <- import("nashpy", convert = FALSE)
# creamos el juego
np1<-np$array(e1)

np2<-np$array(e2)

jg<- nash$Game(np1,np2)
#obtenemos el equilibrio 
equilibria <- jg$support_enumeration(non_degenerate = FALSE)
for (i in iterate(equilibria)){

}
#el último i contiene la estrategia mixta, si la hay
  return(list('p' = py_to_r(i)[[1]][1],
              'q' = py_to_r(i)[[2]][1]))
}

# esta función resuelve el problema completo
# devolviendonos la solución del mismo 

solucion <- function(k = 400,
                     m = 2000,
                     l_ps = 0.8,
                     l_xb = 0.9) {
  # calculamos los diferentes ejercicios de Bertrand
  precios1=matrix(
        nrow = 2,
        ncol = 2,
        dimnames = list(
          c("No hace crossplay", "Hace crossplay"),
          c("No hace crossplay", "Hace crossplay"))
        )
  precios2=precios1
  beneficios1=precios1
  beneficios2=precios1
  for (i in 0:1) {
    for (j in 0:1){
      r= resolucion(
        k = k,
        m = m,
        c_ps = j ,
        c_xb = i,
        l_ps = l_ps,
        l_xb = l_xb
      )
    precios1[j+1,i+1]=r$precios[1]
    precios2[j+1,i+1]=r$precios[2]
    beneficios1[j+1,i+1]=round(r$beneficios[1])
    beneficios2[j+1,i+1]=round(r$beneficios[2])
    }}
  # usamos el programa anterior para calcular 
  # la solución del problema 2x2
  solution2x2 = solucion2(precios1, precios2)
  p = solution2x2$p
  q = solution2x2$q
  # nos devuelve la solución junto con los 
  # beneficios esperados
  return(
    list(
      'p' = p,
      'q' = q,
      'beneficiosesperados' = list(
        'ps' = cbind(p, 1 - p) %*%
          beneficios1 %*%
          c(q, 1 - q),
        'xb' = cbind(p, 1 - p) %*%
          beneficios2 %*%
          c(q, 1 - q)
      ),
      'beneficiosps' = as.data.frame(beneficios1),
      'beneficiosxb' = as.data.frame(beneficios2),
      #obtenemos los precios que pone cada jugador en cada opción
      'preciosps' = as.data.frame(precios1),
      'preciosxb' = as.data.frame(precios2)
    )
  )
}

sol =solucion()

sol$beneficiosesperados

## $ps
##            [,1]
## [1,] 3665411948
## 
## $xb
##            [,1]
## [1,] 4376168454

sol$p

## [1] 0

sol$q

## [1] 0

sol$preciosps

##                   No hace crossplay Hace crossplay
## No hace crossplay          420.3572       249.2645
## Hace crossplay             440.2721       439.8871

sol$preciosxb

##                   No hace crossplay Hace crossplay
## No hace crossplay          408.0221       413.5474
## Hace crossplay             299.0318       460.4085

sol$beneficiosps

##                   No hace crossplay Hace crossplay
## No hace crossplay        3063835537      460007800
## Hace crossplay           3605449303     3665411948

sol$beneficiosxb

##                   No hace crossplay Hace crossplay
## No hace crossplay        2747973569     2749898042
## Hace crossplay            916862215     4376168454

ANÁLISIS DE SENSIBILIDAD

# Función que realiza un mallado 
# de la k y la m 
sensitivity <- function(lstk, lstm)
{
    #usamos las variables auxiliares 
    psolaux = c()# almacena las p obtenidas
    qsolaux = c()# almacena las q obtenidas
    pssolaux = c()# almacena los beneficios de 
    # ps obtenidos
    xbsolaux = c()# almacena los beneficios de 
    # xb obtenidos
    ps2solaux = c()# almacena los precios de 
    # ps obtenidos para c_ps=0 y c_xb=1
    xb2solaux = c()# almacena los precios de 
    # xb obtenidos para c_ps=0 y c_xb=1
    for (i in lstm)# bucle para la m
    {
        for (j in lstk)# bucle para la k
        {   # obtenemos la solución
            sol = solucion(k = j, m = i)
            # almacenamos los datos que interesan
            psolaux = c(psolaux, round(sol$p,
                4))
            qsolaux = c(qsolaux, round(sol$q,
                4))
            pssolaux = c(pssolaux,
                round(sol$beneficiosesperados$ps,
                  4))
            xbsolaux = c(xbsolaux,
                round(sol$beneficiosesperados$xb,
                  4))
            ps2solaux = c(ps2solaux,
                round(sol$preciosps[2,
                  2], 2))
            xb2solaux = c(xb2solaux,
                round(sol$preciosxb[2,
                  2], 2))
        }
    }
    # transformamos los datos a dataframes
    plista = as.data.frame(matrix(psolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    qlista = as.data.frame(matrix(qsolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    pslista = as.data.frame(matrix(pssolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    xblista = as.data.frame(matrix(xbsolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    xb2lista = as.data.frame(matrix(xb2solaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    ps2lista = as.data.frame(matrix(ps2solaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    # devolvemos los diferentes dataframes
    return(list('p' = plista, 'q' = qlista,
        'beneficiosps' = pslista, 'beneficiosxb' = xblista,
        'preciosps' = ps2lista, 'preciosxb' = xb2lista))
}

# creamos el mallado
lstk = seq(200, 600, length.out = 9)
lstm = seq(1000, 3000, length.out = 9)
# evaluamos el mallado 
lista = sensitivity(lstk, lstm)

# obtenemos todas las tablas

tabla(lista[["p"]])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	0	0	0	0	0.8752	0.848	0.8250	0.8056	0.7898
k= 250	0	0	0	0	0.0000	0.000	0.8726	0.8519	0.8328
k= 300	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.8721
k= 350	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 400	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 450	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 500	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 550	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 600	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000

tabla(lista[["q"]])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	0	0	0	0	0.9928	0.977	0.9651	0.9548	0.9470
k= 250	0	0	0	0	0.0000	0.000	0.9920	0.9796	0.9689
k= 300	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.9917
k= 350	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 400	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 450	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 500	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 550	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000
k= 600	0	0	0	0	0.0000	0.000	0.0000	0.0000	0.0000

tabla(lista[['beneficiosps']])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	1507584083	1965591008	2380143772	2752750928	2652873864	2906595062	3140921257	3354637132	3551799715
k= 250	1546131330	2041127306	2501788127	2926159824	3314630443	3669686131	3444385258	3695436686	3929660465
k= 300	1568107579	2087385417	2579991822	3041689354	3472337098	3871936572	4241712304	4585411821	4236637708
k= 350	1581317042	2116865447	2631766986	3120899396	3583849435	4018537734	4426469866	4808418765	5165654640
k= 400	1589610198	2136159879	2667660319	3178076864	3665411948	4127766690	4566005972	4979761959	5370182729
k= 450	1594821540	2149494742	2692764130	3218917373	3725420728	4210351901	4672534439	5113045012	5532264982
k= 500	1598154405	2158630791	2710735907	3249232976	3770762571	4273313057	4755872762	5218676287	5660733657
k= 550	1600444871	2164962963	2724059154	3271892288	3805303363	4322349683	4820988533	5302190893	5765104296
k= 600	1601867661	2169741342	2733853579	3289212343	3832239148	4360797813	4874107560	5369454961	5849185658

tabla(lista[['beneficiosxb']])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	1839819353	2361427842	2830308934	3249545230	2194850974	2395607578	2577488075	2744601586	2898366183
k= 250	1901269854	2468318160	2991746669	3472061795	3909612095	4308629825	2858576629	3059381206	3243179154
k= 300	1941150950	2539320374	3101782734	3626632952	4112964205	4564244623	4979848153	5365708599	3527160129
k= 350	1968901134	2588896408	3179063504	3737922487	4262993529	4754425596	5215309730	5645356619	6046917677
k= 400	1989055132	2625235912	3237404138	3821430326	4376168454	4900972430	5397094732	5864128281	6304138902
k= 450	2004316972	2652984610	3281296274	3885561098	4464085862	5015532028	5540123887	6039106288	6512146528
k= 500	2016297824	2674562585	3315837431	3936115123	4533714984	5106978440	5655492205	6180745235	6681047097
k= 550	2025885318	2691853692	3343519388	3976700336	4589626178	5181138827	5749669826	6296532312	6821040583
k= 600	2033627763	2706077134	3366207940	4010024095	4635529465	5242055425	5827721737	6392924016	6938014217

tabla(lista[['preciosps']])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	262.54	290.32	312.87	331.48	346.95	359.97	371.05	380.59	388.87
k= 250	274.02	306.94	334.62	358.20	378.23	395.51	410.39	423.46	434.87
k= 300	282.32	319.34	351.29	379.04	403.09	424.29	442.75	459.27	473.89
k= 350	288.62	328.83	364.15	395.57	423.31	447.84	469.87	489.49	507.09
k= 400	293.51	336.33	374.72	409.05	439.89	467.57	492.61	515.12	535.52
k= 450	297.44	342.45	383.21	420.17	453.75	484.18	511.85	537.17	560.23
k= 500	300.68	347.47	390.30	429.49	465.45	498.32	528.41	556.22	581.63
k= 550	303.40	351.69	396.29	437.39	475.37	510.47	542.75	572.72	600.44
k= 600	305.64	355.32	401.43	444.19	483.92	520.97	555.28	587.21	617.02

tabla(lista[['preciosxb']])

	m= 1000	m= 1250	m= 1500	m= 1750	m= 2000	m= 2250	m= 2500	m= 2750	m= 3000
k= 200	273.54	300.44	322.04	339.78	354.46	366.76	377.25	386.25	394.13
k= 250	287.86	320.31	347.22	369.85	389.01	405.43	419.60	431.91	442.73
k= 300	298.52	335.58	367.05	393.95	417.26	437.50	455.14	470.86	484.71
k= 350	306.74	347.59	382.89	413.54	440.64	464.39	485.49	504.27	521.04
k= 400	313.27	357.20	395.95	430.06	460.41	487.34	511.54	533.22	552.76
k= 450	318.52	365.17	406.73	443.80	477.12	507.07	533.99	558.45	580.76
k= 500	322.81	371.77	415.78	455.54	491.49	524.07	553.66	580.73	605.37
k= 550	326.50	377.29	423.56	465.59	503.91	538.93	570.82	600.27	627.40
k= 600	329.59	382.16	430.22	474.35	514.80	551.95	586.29	617.64	646.92

# Función que realiza un mallado
# de la l_ps y la l_xb
sensitivity2 <- function(n,a,b)
{
  lstps=rbeta(n,a,b)
      #usamos las variables auxiliares
    psolaux = c()# almacena las p obtenidas
    qsolaux = c()# almacena las q obtenidas
    pssolaux = c()# almacena los beneficios de
    # ps obtenidos
    xbsolaux = c()# almacena los beneficios de
    # xb obtenidos
    ps2solaux = c()# almacena los precios de
    # ps obtenidos para c_ps=0 y c_xb=1
    xb2solaux = c()# almacena los precios de
    # xb obtenidos para c_ps=0 y c_xb=1
    for (l_ps in lstps)# bucle para la l_ps
        {   # obtenemos la solución
            sol = solucion(l_ps = l_ps)
            # almacenamos los datos que interesan
            psolaux = c(psolaux, round(sol$p,
                4))
            qsolaux = c(qsolaux, round(sol$q,
                4))
            pssolaux = c(pssolaux,
                round(sol$beneficiosesperados$ps,
                  4))
            xbsolaux = c(xbsolaux,
                round(sol$beneficiosesperados$xb,
                  4))
            ps2solaux = c(ps2solaux,
                round(sol$preciosps[2,
                  2], 2))
            xb2solaux = c(xb2solaux,
                round(sol$preciosxb[2,
                  2], 2))
        }

    
    # devolvemos los diferentes dataframes
    return(list('p' = psolaux, 'q' = qsolaux,
        'beneficiosps' = pssolaux, 'beneficiosxb' = xbsolaux,
        'preciosps' = ps2solaux, 'preciosxb' = xb2solaux))
}

set.seed(2)
n=100
a=8
b=2
lista2 = sensitivity2(n,a,b)

#obtenemos todas las tablas

mean(lista2[["p"]])

## [1] 0

mean(lista2[["q"]])

## [1] 0

histo<-function(lst){
  m=mean(lst)
  mo=mlv(lst,method = 'shorth')
  
  q1=quantile(lst, 0.025)
  q2=quantile(lst, 0.975)
  
  hist(lst,freq=FALSE)
  abline(v =m,col = "blue",lwd = 2)
  abline(v =q1,col = "red",lwd = 2)
  abline(v =q2,col = "red",lwd = 2)
  abline(v =mo,col = "green",lwd = 2)
  return(paste(
    'El cuantil para 0.025 es ' , as.character(q1),'
    El cuantil para 0.975 es ',as.character(q2),'
    La media es ' ,as.character(m),'
    La moda es ' ,as.character(mo),'
    La desviación típica es ',as.character(sd(lst))
  ))
}

# analizamos los beneficios de Playstation
histo(lista2[['beneficiosps']])

## [1] "El cuantil para 0.025 es  1615273846.525 \n    El cuantil para 0.975 es  5198874991.925 \n    La media es  3802338582.66 \n    La moda es  4302440780.54 \n    La desviación típica es  1006283136.00399"

# analizamos los beneficios de Xbox
histo(lista2[['beneficiosxb']])

## [1] "El cuantil para 0.025 es  3883555405.2 \n    El cuantil para 0.975 es  4605210455.425 \n    La media es  4374980663.61 \n    La moda es  4482525143.96 \n    La desviación típica es  206616743.334366"

histo(lista2[['preciosps']])

## [1] "El cuantil para 0.025 es  350.83675 \n    El cuantil para 0.975 es  480.531 \n    La media es  439.4743 \n    La moda es  458.7956 \n    La desviación típica es  37.2257921408686"

histo(lista2[['preciosxb']])

## [1] "El cuantil para 0.025 es  447.37075 \n    El cuantil para 0.975 es  466.60975 \n    La media es  460.3759 \n    La moda es  463.4226 \n    La desviación típica es  5.48209631344607"